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I. INTRODUCTION 

A. Background 

Recently Soltan-Panahi et al. found evidence of 
a zero quasi-momentum "Twisted Superfiuid (TSF) 
phase" of a two-component Bose-Einstein conden- 
sate (BEC) trapped in a spin-dependent honeycomb 
lattice 1 . A TSF is characterized by Bose-Einstein 
condensation into a state whose order parameter (a 
macroscopically occupied single particle wavefunction) 
has a spatially varying phase. The simplest example is 
condensation at finite momentum. Alternatively, in a 
non-Bravais lattice where the unit cell involves multiple 
sites, one can have a TSF at zero quasi-momentum if 
the phase of the order parameter varies throughout the 
unit cell. We model Soltan-Panahi et al.'s experiment 
[T] with a mean field Gross-Pitaevskii function. We 
find that the TSF phase is absent within mean field 
theory thus suggesting that the observations are due to 
non-mean field effects. 

Twisted Superfiuid phases are quite exotic; the 
phase twists are naturally associated with microscopic 
currents. Moreover, the present example involves spon- 
taneous symmetry breaking, and provides a setting for 
studying phase transition physics. Analogous physics 
can be found in magnetic systems [2] and in the excited 
states of lattice bosons [3J 2] . 

To see the connection between TSF and currents, con- 
sider a hopping model on a honeycomb lattice : 

H = -t^alaj +ajai, (1) 

<ij> 

where a; annihilates a boson on site i, t is the en- 
ergy scale of hopping and (ij) represent nearest neighbor 
sites. The rate of change of the density on site i is : 

r i = ^(n i )= J2 -t Im < Sj&j >, (2) 

< neighbours j> 
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and one interprets — t Im < a[a,j > as the net current 
flowing from site j to site i. In the case of a condensed 
system, the operators aj can be replaced by c-numbers 
and any non-trivial phases give rise to currents. 

The experimental evidence suggests that the order pa- 
rameter is uniform on each of the triangular sub-lattices 
of the honeycomb lattice : 

(a j± ) = Vh^cxp(±i 5/2), (3) 

where the + and - signs are taken on the two dis- 
tinct sub-lattices. For the simple model in eq.(l), 
j^rii — ±3 t ^/n + n_ Sin((5), and the density is not 
stationary unless 5 = or w. This argument misses 
the role of interactions and is therefore incomplete. In 
particular, interactions can give a Sin(2<5) contribution 
to Ti, and yield steady state for (5 ^ or it. This can 
be interpreted as counter-flowing single particle and 
pair currents which yield zero net flow (see Section V). 
Regardless of the interpretation, the sub-lattice phase, 
5 is worthy of study. 



B. Experimental Evidence 

Soltan-Panahi et aVs evidence for non-trivial phases 
comes from time-of-flight (TOF) expansion, a technique 
where all trapping fields are removed and the atomic en- 
semble falls freely under gravity. Neglecting interactions 
[7], the long-time real space density profile is simply 
the initial density in momentum space. For the special 
case of a BEC, n k = \ip(k)\ 2 = \J cxp(+ik.r)^(r)| 2 
and the time-of-flight density profile is the squared 
modulus of the Fourier transform of the order param- 
eter, ^(r). As schematically illustrated in Figure 1, if 
ip(r) is real, and has the symmetry of the honeycomb 
lattice, its Fourier transform (and consequently the 
TOF pattern) is six fold symmetric. This symmetry 
persists even if the densities on the two sub-lattices 
differ, forming a three-fold symmetric charge density 
wave (CDW) as illustrated in Figure 2. This six- fold 
rotational symmetry of the TOF pattern is a conse- 
quence the point group symmetry of the lattice (C^y) 
and the relation k) = V*(k), which holds for 
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FIG. 1. Schematic of the Time-of-Flight (TOF) pattern for 
a superfluid in a 2D honeycomb lattice. Larger darker dots 
correspond to more particles with a given momentum. The 
complex numbers t and z represent the amplitudes of the 
Fourier transform of the condensate wavefunction at k — 
(f,0) and k = (^r,^)The twisted superfluid (TSF) is 
described by jtj 7^ |z|. 



real ^(r). Therefore, if the TOF pattern becomes 
three-fold symmetric (|t| 7^ |z| in Figure 1), the real 
space wavefunction has different phases on the two sub- 
lattices. The experimentalists see exactly this signature. 

Soltan-Panahi et al. work with a two-component mix- 
ture of 87 Rb atoms in a spin-dependent honeycomb lat- 
tice PQ. The two spin states can form out-of-phase 
CDWs. As was emphasized in [T] this density wave does 
not lead to the breakdown of the rotational symmetry 
in time-of- flight. The symmetry breaking was opposite 
for the two species (i.e j|4 = l|4, where t, z denote the 
amplitudes of the TOF peaks in Fig. 1). 



II. THE MODEL 

Within a mean field model, we will investigate the 
relative stability of twisted or ordinary superfiuids. The 
energy of a two component BEC, described by macro- 
scopic wavefunctions ipi and ^2 is : 



E3D 



2m 



|VVv(r)| 2 + Vv(r)|Vv(r)| 5 



77°" 

er = l,2 ^ 



+ V conl (r)(\M*)\ 2 + \Mr)\ 2 ) 



(4) 



Here, U£ D = 4t ^°° is the intra-species interaction 
energy (a^ is the intra-species scattering length for 



FIG. 2. The density wave formed in a honeycomb lattice 
for the rr%F = 1 atoms. The points represent lattice sites. 
Larger points indicate more atoms. This pattern is periodi- 
cally repeated. A complementary density wave is formed by 
tuf = — 1 atoms. This density wave does not lead to a 6-fold 
symmetry breaking in TOF unless additional phases appear 
on the sites. 



species a), while W$d = 4ll ' m '" i is the inter-species 
interaction energy (di2 is the inter-species scattering 
length). We focus on the case in Q], where the states 
1 (described by ipi) and 2 (described by ^2) are the 
|F = l,m F = 1) and |F = l,m F = —1) states of 
87 Rb. For these two hyperfine states of 87 Rb atoms, 
U^ D ,U^ D a.nd W are almost equal (« lOOao where 
ao is the Bohr radius). In principle collisions can 
connect these hyperfine states to others (for example 
|F = l,mp = 0)). For the experimental parameters, 
these processes are off-resonant and the two-component 
Bose gas model describes the physics. 



In the experiment pQ, the honeycomb lattice is gen- 
erated by 3 lasers yielding a potential V^(r) = Vh GX (r) ± 
aB c ff (r) where, state 1 sees the sign '+' and state 2 sees 
the sign '-' (with a — 0.13) and 



B, 



x (r) = 2 Vl at (Cos[k L bi-xj+CosIkL b 2 .x] 
+ Cos[k L b 3 .x]) 

W(r) = 2V3 M at (Sin[k L b!.x]+Sm[k L b 2 .x] 
+ Sin[k L b 3 .x]) 



(5) 



(0) 
(7) 
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where, bl = — ^e x — ^e y ; b2 = e x ; b3 

and kL = 2^/^-k/\i j (\l is the laser wavelength and is 
830 nm for the experiment under discussion). With 
these considerations Vi a t is the height of the barrier 
between neighboring sites. The difference between the 
maximum and minimum values of Vhcx(r) is 4 T4at- 
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The experiment uses a separate set of lasers to provide 
strong confinement in the third dimension, V con f(r): 



III. METHOD 



27T 

V coni (r) = Vid Cos[- — z] 
Aid 



VlD / 27T ■ 2 ^2 

2 Aid 



(8) 



This potential restricts the dynamics to two dimensions 
and we may take the wavefunction of the BEC in the 
third direction to be constant and Gaussian. Then the 
energy can be written as : 



E 2D = f d 2 r h|Vvv(r) + K(rM(r) 

/ z — ' Z/77 

° -r 1 1 



+ ^l^(r)| 4 ]+^ 2C |^(r)| 2 |^(r)| 2 (9) 



where U 2D = lW Vm ^ - and W, 



'2D — 



WsDyJ^fnri^ 21 - In the ex P eriment CO, Aid = Aj, = 
830 nm and V\d — 8.8Er. For these parameters, the 
weakest lattice yielding a Mott state is Vut ~ 3.5Er for 
two particles per unit cell within the Gutzwiller mean 
field approximation [S]. 

From the time-of-flight images obtained, the ex- 
perimentalists observe a breakdown of the six-fold 
rotational symmetry in momentum space for lattice 
depths, Mat ranging from about 1 to 4 Er, where 



E R = 



2m Rb \ 2 L 



(mRb is the mass of Rb atoms) 



We assume a form of - 0i( r ) and V'aCr) which is con- 
sistent with the timc-of-flight measurements : 

V>i(r) =^Vi(k)exp(-i k.r), (10) 

k 

^(r) = ^ 2 (k)exp(-ik.r). (11) 

k 

where k are the reciprocal lattice vectors of a honey- 
comb lattice. We insert this variational ansatz into 
eq.Q and minimize the energy with respect to the 
set of variational parameters ^i(k) and 02 (k)- We 
find from our simulations that for all experimental 
parameters ipi(k) = -02 (k), where 0>|(k) is the complex 
conjugate of 02 (k). This result is sensible and implies 
ipi and ip2 are related by a lattice translation. 

We perform the variational minimization in Fourier 
space rather than real space (where such minimiza- 
tion is usually done). This is equivalent to solving 
the Gross-Pitaevskii equation in real space within a 
single unit cell with periodic boundary conditions. 
Computationally, we find momentum space to be more 
efficient. Moreover, the experimental probes are all in 
momentum space. 



In k-space, the energy expression, eq.(|9]) becomes 



E 2 d 
Er 



_Et 3 *V«*(k)Vi(k) 

{k,ki,k 2 ,k 3 }e£ l=1 . 2 

M(k 1 )^(k 2 )^(k 2 -k 1 ) 

|vr(ki)^(k2)^(k 3 )^(ki + k 2 - 
WVi*(ki)Vi(k 2 )V'2(k3)V'2(ki + k 3 



ks)] 



(12) 



where C stands for the reciprocal lattice i.e 
k = (ajbi + a2b2), where ai and &2 are inte- 
gers. One can also generate this lattice from one of 
bl, b2 and b3, all explicitly given following Eq.(7). 
All energies (Vi,U and W) are adimensionalized by Er. 

While we carried out unrestricted minimizations, our 
results are best illustrated by considering an ansatz 
where the low momentum physics is characterized by 
2 complex numbers t and z. In particular, we take 
V»i(k) = t and -0 2 (k) = z for k = {bl, b2, b3} and 
V>i(k) = z and ip 2 (k) = t for k= {-bl, -b2, -b3}. 
In terms of their real and imaginary parts, we write 



t = t r + i tj and 

Z = Z r + i Zj. 



(13) 
(14) 



The order parameter (OP) for the TSF phase is given 
by: 



OP = 



|t| 2 



M 2 + t 



(15) 



By construction, OP has a non-zero value in the 
TSF phase and is zero for a uniform condensate. 
Soltan-Panahi et al. measure this order parameter in 
their experiment to demonstrate the presence of the 
TSF phase. 

For our minimization, we restrict ourselves to |k| < 6 
giving us 159 complex variational parameters. We find 
that there are no differences if we use |k| < 4 instead 
(62 variational parameters). Therefore, we believe our 
results faithfully reflect what would be found if an 
infinite number of Brillouin zones were included. We 
gain further confidence in the convergence of our results 
by noting that the fraction of population occupying 
the |k| = 5 state when U = 0.05£ fl and V\ &t = 3.8E R 
is about 0.1%. In our simulations, we vary U in the 
range 0.04_Er to U — 0.2Er corresponding to various 
strengths of the transverse confinement. We also vary 
a in the range 0.08 to 0.3, corresponding to varying 
amounts of detuning of the laser beams. 
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IV. RESULTS 

We do not find any evidence for the existence of the 
Twisted Superfluid phase despite an extensive search 
of the parameter space. Since Eq.(12| is a quartic 



form, it will in general have multiple minima and a 
number of other stationary points. The most grave 
concern with our results is that we might not have 
found the global minimum. To some extent, we can 
alleviate this concern by noting that the experiment 
finds a continuous symmetry breaking as a function of 
lattice depth. It therefore suffices to establish that our 
solution is a dynamically stable local minimum which 
is continuously connected to the symmetry-unbroken 
ground state at V\ a t = 0. 
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A. Local Energetic Stability 

We check whether whether we have found a true min- 
imum by looking at the eigenvalues of the Hessian. The 
Hessian H is a matrix whose elements are given by : 

d 2 E 
ddiddj ' 

where aj and Oj are real variational parameters (cor- 
responding to the real and imaginary parts of ?A(k)). 
We find that for all parameters, the eigenvalues of 
H axe positive. This implies that we have at least 
found a local minimum. In Figure [3j we plot the 
minimum eigenvalues of the Hessian for different values 
of the lattice depth (Viat) at the illustrative interaction 
strength, U — 0.05Er and a — 0.14, for five particles 
(of each species) per unit cell. 

We further illustrate the stability of our theory by 
doing two separate numerical experiments : 

(a) Fix the ratio of z r (Re[z]) to t r (Re[t]) and vary 
the remaining variational parameters to find the energy 
minima. We find that the minimum of the energy 
occurs when z r : t r — 1 and there are no other local 
minima. The dotted curve shows this in Figure |4j 

(b) Fix the ratio of z; (Im[z]) to tj (Im[t]) and vary 
the remaining variational parameters to find the energy 
minima. We find that the minimum of the energy 
occurs when %- x : tj =1 and there are no other local 
minima. The solid curve shows this in Figure [4j 

We conclude that there is no second order phase tran- 
sition within mean field theory. 



B. Local Dynamic Stability 

We also check whether the minima found is unstable 
against perturbations. This is done by looking at the 



FIG. 3. Minimum eigenvalue of the Hessian, Ao in the Nor- 
mal superfluid phase plotted against the lattice depth, Vi B t 
(in units of Er) when U = 0.05Er and 5 particles (of each 
species) are present per unit cell. All the eigenvalues of the 
Hessian are positive, thereby showing the stability of the nor- 
mal phase. We conclude that there is no Twisted superfluid 
phase for these potential depths. This result is illustrative 
of all parameter ranges we explored. 



Gross-Pitaevskii equation 



dt dip* 



This would imply 



dSaj SE 



Of 



^ ddjdai 



(17) 



(18) 



Taking the real and imaginary parts of both sides, we 
get the eigenvalue equations 



fko u = Mu 



(19) 



where, 



M 



Re[H] 
lm[H] 



-lm[H) 
Re [if] 



We look at the eigenvalues of this matrix, M. A 
complex eigenvalue would signify the presence of a 
mode which will grow with time, thus rendering this 
ground state unstable. We find that all the eigenvalues 
are real. Thus, the minimum that we have found is 
also dynamically stable. This is a generic feature of 
quantum systems: Energetic stability implies dynamic 
stability [5J. 



V. BEYOND MEAN FIELD THEORY 

In order to explain the experimental observations, one 
clearly must go beyond our model. One possible direc- 
tion is including fluctuation effects. These become large 



5 



Energy 
Er 

-20.733 r 



-20.734 
-20.735 
-20.736 
-20.737 
-20.738 
-20.739 
-20.740 



and 




0.8 



0.9 



1.0 



1.1 



1.2 



z:t 



FIG. 4. Slice through the energy landscape at Vj a t = I-SEr 
and U — 0.05-Bij and 5 particles (of each species) per unit 
cell. Dotted curve: The ratio Re [z]: Re [t] is varied and the en- 
ergy is found by minimizing with respect to the other vari- 
ational parameters. Solid curve: Same, but with varying 
Im[z]:Im[t]. We find that the overall energy minimum oc- 
curs when Re[z] = Re[t] and Im[z] = Im[t]. 



as one approaches the Mott transition (which as pre- 
viously explained may occur at lattice depths as low 
as Via.% ~ 3.5-Er for two particles per unit cell). In this 
section, we explore the structure of a more general treat- 
ment and evaluate how a twisted supcrfiuid (TSF) could 
emerge. The most direct approach is is to reduce Eq.(9) 
to a tight-binding model by writing 



n), 



(20) 



where d a i annihilates a boson of species a on site i, at 
position r;, where (j>(r — ri) is the Wannier state wave- 
function. With this substitution, the Hamiltonian be- 
comes : 



full 



H 1 +H 2 + .. 



(21) 



where Hq involves only the onsite terms, H\ involves the 
nearest neighbors, H2 involves the next nearest neigh- 
bors and so on. Our mean field theory includes all of 
the terms in Eq.(|2~ij). At large lattice depths, the Wan- 
nier states are well localized, and each subsequent term 
in Eq.(21) is smaller. In addition to a large number 
of other expressions, Hi contains single particle hop- 
ping and pair hopping terms. Neglecting the spin index, 
these terms are 



H ( sp h) = _ t (gtg. + gtg.) 



ff(P h ) = t'(af at ajaj + 



clj clj ctjSLi 



(23) 



respectively. Within mean field theory, these terms give 
a contribution to the energy 

jff (-ph) +jff (ph) = 3 (-t^EICos(S) + t'n+n_Cos(2(5)) , 
_ (24) 
v / h±exp(±i 8/2), where the + and - 



where (a,j±) 

signs are taken on the distinct triangular sub-lattices. 
This approximation neglects much of the physics in- 
cluded in sections III and IV, but nicely explains a mech- 
anism for developing a TSF. Namely, if (4t'n + n_ > 
t^/n+n-), Eq.(24) is minimized by making 6 ^ (cf. 
Methods section of [1] ) . As illustrated by the absence of 
a TSF in the full mean field theory, the neglected terms 
must suppress this symmetry breaking. One could how- 
ever imagine that Mott physics may renormalize param- 
eters and change this result. For example, Mott physics 
will send ij( s P h ) to zero in the insulating phase, while 
having a smaller effect on terms such as countcrflow 
hopping: 



T a j; a jT a i; 



a i4 a jt^j4%) 



(25) 



<ij> 



(22) 



If such terms are responsible for the TSF, then the 
experimental observations can be interpreted as a 
signature of "transverse" magnetic ordering or at least 
short range magnetic correlations. 

One could experimentally explore the importance 
of quantum fluctuations by changing the strength of 
the transverse confinement, V con f(r). This impacts the 
onsite interaction while leaving the hopping largely 
unaffected. Systematically studying the order pa- 
rameter as a function of lattice depth and transverse 
confinement may provide the key to understand this 
phenomenon. 

Another direction worth exploring is the modeling of 
the time-of-flight expansion or the turn-off of the optical 
lattice. Phases accumulated during the expansion could 
mimic a TSF. Estimates of these phases suggest however 
that they are too small [7]. 



ACKNOWLEDGEMENTS 

This paper is based on work supported by the Na- 
tional Science Foundation under Grant no. PHY- 
1068165. 



[1] P. Soltan-Panahi el at, Nature Physics 8, 71-75 (2012). [3] G. Wirth, M. Olschlager and A. Hemmerich, Nature 
[2] J. Struck et ai, Science 333, 996 (2011). Physics 7 ,147-153 (2011). 



G 



[4] M. Olschlager, G. Wirth and A. Hemmerich Phys. Rev. 

Lett. 106, 015302 (2011). 
[5] P. Soltan-Panahi et al, Nature Physics 7, 434-440 (2011). 



[6] C. J. Pethick and H. Smith, Bose-Einstein Condensation 
in Dilute Gases. Cambridge University Press, Cambridge 
(2002). 

[7] J. N. Kupferschmidt and E. J. Mueller , Phys. Rev. A 
82, 023618 (2010). 



